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For the recently introduced algorithms to solve the time-dependent Maxwell equations Jy], we 
■ construct a variable grid implementation and an improved spatial discretization implementation that 

Q preserve the exceptional property of the algorithms to be unconditionally stable by construction. 
We find that the performance and accuracy of the corresponding algorithms are significant and 
illustrate their practical relevance by simulating various physical model systems. 
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Oh- I. INTRODUCTION 

In a recent paper, we introduced a family of algorithms to solve the time-dependent Maxwell equations (TDME) M . 
Salient features of these algorithms include the rigorously provable unconditional stability for e?-dimensional systems 
(d = 1,2,3) with spatially varying permittivity and permeability, as well as the exact conservation of the energy density 
of the electromagnetic (EM) fields. Furthermore, we have demonstrated that - without affecting the unconditional 
O stability of the algorithms - the order of accuracy in the time integration can be systematically increased. An important 
aspect that has not been considered in our earlier work concerns the effect of the discretization of space on the 
accuracy of the algorithms. Previously, we employed only the most simple spatial discretization, namely the central- 
4^ ' difference scheme on a cartesian grid with a constant mesh size jjj . We refer to this basic spatial discretization scheme 
as the poor man's implementation. Many numerical simulations of realistic physical systems require algorithms with a 
more accurate spatial discretization and a more flexible spatial grid for an optimal use of computer resources (CPU time 
t-H , and computer memory) . In the present paper we show that implementing a fourth-order accurate approximation of the 
spatial derivatives and a spatial grid of variable mesh sizes preserve the unconditional stability of the algorithms. We 
simulate various physical model systems using these new implementations to demonstrate the significant improvement 
with respect to the required computer resources in the computation of eigenmode spectra and to study systematically 
' the temporal and spatial accuracy of the algorithms. 

Our presentation is organized as follows: We recapitulate the theory of constructing unconditionally stable algo- 



rithms to solve the TDME in Sec. II and describe the basic properties of the poor man's implementation in Sec. IIP 
Then, in Sec. IV and Sec. [V| we present the implementation of, respectively, the variable grid and the improved spatial 



^5 ■ discretization. Our conclusions are given in Sec. VI 

o : 
•1— 1 

^ II. UNCONDITIONALLY STABLE ALGORITHMS TO SOLVE MAXWELL'S EQUATIONS 

^ : 

We consider a d-dimensional model system of EM fields in a medium with spatially varying permittivity and/or 
. permeability, surrounded by a perfectly conducting box. In the absence of free charges and currents, the EM fields in 
such a system satisfy Maxwell's equations 

— H = VxE and — E=-VxH, (1) 

at \i at e 

diveE = and divH = , (2) 
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where H = (H x (r, t), H y (r, i), H z (r, t)) T and E = (E x (r, t), E v (r, t), E z (r, t)) T denote, respectively, the magnetic and 
the electric field vector. The permeability and the permittivity are given by /i = /i(r) and e = e(r). For simplicity of 
notation, we will omit the spatial dependence onr= (a;, y, z) T unless this leads to ambiguities. On the surface of the 
perfectly conducting box the EM fields satisfy the boundary conditions j| 

n x E = and n ■ H = , (3) 

with n denoting the vector normal to a boundary of the surface. The conditions Eqs. (j^) assure that the normal 
component of the magnetic field and the tangential components of the electric field vanish at the boundary Q . Some 
important symmetries of the Maxwell equations ([j])-(Q) can be made explicit by introducing the fields 

X(t) = s/JM{t) and Y(t) = VeE(t) . (4) 

In terms of the fields X(t) and Y(t), the TDME (g) read 

d_(X(t)\_( _ -^Vx^/ X (i)^^X(i) 



di\Y(t)J-y^ X ^ ■ )\Y(t)J= H {Y(t)l- (5) 
Writing *(t) = (X(t), Y(t)) T , Eq. (§) becomes 

^*(i)=W*(t). (6) 

It is easy to show that 7i is skew-symmetric, i.e. Ti. T = —H, with respect to the inner product (vP|\|j') = J v \I/ T • \&' dr, 
where V denotes the volume of the enclosing box. The formal solution of Eq. (jfy is given by 

= E7(t)¥(0) = e m *(0), (7) 

where Vt'(O) represents the initial state of the EM fields. The operator [/(t) = e* w determines the time evolu- 
tion. By construction ||*(i)l| 2 = (*(*)!*(*)) = Jy [eE 2 (t) + ^H 2 (i)] dr, relating the length of *(t) to the en- 
ergy density w(t) = eE 2 (t) + fiU 2 {t) of the EM fields @. As U{tf = U(-t) = U' 1 ^) = e~ m it follows that 
(f/(f)*(0)|f7(t)*(0)) = (*(t)|*(i)) = (*(0)|*(0)). Hence the time-evolution operator ?7(t) is an orthogonal trans- 
formation, rotating the vector *&(t) without changing its length ||^||. In physical terms this means that the energy 
density of the EM fields does not change with time, as expected on physical grounds 0. 

A numerical procedure that solves the TDME necessarily starts by discretizing the spatial derivatives. This maps 
the continuum problem described by TL onto a lattice problem defined by a matrix H . Ideally, this mapping should 
not change the basic symmetries of the original problem. The underlying symmetry of the TDME suggests to use 
matrices H that are real and skew-symmetric. Since formally the time evolution of the EM fields on the lattice 
is given by *&(t + r) = J7(r) 1 4'(t) = e r - ff \I/(t), the second ingredient of the numerical procedure is to choose an 
approximation of the time-evolution operator U(t). The fact that U(t) is an orthogonal transformation is essential 
for the development of an unconditionally stable algorithm to solve the Maxwell equations [|J . A systematic approach 
to construct orthogonal approximations to matrix exponentials is to make use of the Lie- Trotter-Suzuki formula [|[ |^] 

e t(H 1 + ...+H p ) = Hm f~TT e tHi/m\ ^ (g) 
\i=l / 

and generalizations thereof || |) . Applied to the case of interest here, the success of this approach relies on the basic 
but rather trivial premise that the matrix H can be written as H = 5^<=i Hit where each of the matrices Hi is real 
and skew-symmetric. Expression Eq. suggests that 

U 1 (T) = e THl ...e TH - (9) 

might be a good approximation to U(t) if r is sufficiently small. In fact, it can be shown that U(t) and U\{t) are 
identical up to first order in r 0. Most importantly, if all the Hi are real and skew-symmetric, U\{t) is orthogonal 
by construction. Therefore, by construction, a numerical scheme based on Eq. (^[) will be unconditionally stable. The 
product-formula approach provides simple, systematic procedures to improve the accuracy of the approximation to 
U(t) without changing its fundamental symmetries. For example the orthogonal matrix 

U 2 (r) = U?(-t/2)Ui(t/2) (10) 
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is identical to U(t) up to second order in r p, ||. Suzuki's fractal decomposition approach H gives a general method 
to construct higher-order approximations based on U\{t) or U 2 (t). A particularly useful approximation, which is 
identical to U(t) up to fourth order in r, is given by H 

U 4 (t) = U 2 {aT)U 2 {aT)U 2 ((l - Aa)T)U 2 {aT)U 2 {aT) , (11) 

where a — 1/(4 — 4 1 / 3 ). From Eqs. (p|)-(pd|) it follows that, in practice, an efficient implementation of a scheme based 
on Ui(t) is all that is needed to construct the higher-order algorithms Eqs. ( |l0|) and (f[l]). In many applications the 
approximations U u (t) to the time-evolution operator U(t) have proven to be very useful j|, [| [?], || ||, [l(], [II], [1^, [D], 
|lj, |l5|] and turn out to be equally useful for solving the TDME Q. In particular, it can be shown that the difference 
between the exact EM field vector \I>(i) = ?7(i) 1 3>(0) and the approximate one, *„(<) = f7„(i)vP(0), is bounded by 

\\(U(t)-U n (t))*(0)\\ = ||*(i)-* B (t)|| < C n tr n , (12) 

where C n is a constant. The rigorous upper bound on the error of the EM field vector will be used to specify 
unconditionally stable algorithms by the temporal and spatial accuracy of the computed EM field. We denote an 
algorithm by TnSm if its implementation involves a time integration based on U n {r) and a spatial discretization 
scheme based on an mth-order accurate approximation of the spatial derivatives. 



III. POOR MAN'S IMPLEMENTATION 



In this section, we briefly recapitulate the construction of the unconditionally stable algorithm to solve Maxwell's 
equations in a one-dimensional (ID) system. Furthermore, we discuss general properties of this implementation 
refering also to the two-dimensional (2D) and three-dimensional (3D) case. 

Maxwell's equations for a ID system extending along the a;-axis contain no partial derivatives with respect to y 
or z. Also e and /i do not depend on y or z. Under these conditions, the TDME reduce to two independent sets of 
first-order differential equations The solutions to these sets are known as the transverse electric (TE) mode and 
the transverse mag netic (TM) mode j|. Restricting our considerations to the TM-mode, it follows from Eq. (||) that 
the magnetic field H y {x, t) = X y (x, t)/y/ fi(x) and the electric field E z (x, t) = Y z (x, t)/y/e(x) are solutions of 



0_ 



Xy(x, t) 



1 



d Y z {x,t) 



1 d ( X y {x,t) 



Note that in ID the divergence of H y (x,t) and E z (x,t) is zero, hence Eqs. 
central-difference scheme, which yields a second-order accurate approximation of the spatial derivatives, we obtain 

1 fY z {i+l,t) Y z {i-l,t) 



(13) 
(14) 

are automatically satisfied. Using the 



Wt Y ^ 
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Xy(j+l,t) Xy(j-l,t) 



(15) 
(16) 



where the spatial coordinate of an EM field component is specified through the lattice index i, e.g . X y (i,t) stands 
for X y (x — iS/2,t), and 5/2 the distance between two neighboring lattice points. Following Yee [|16| it is convenient 
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FIG. 1: Positions of the two TM-mode EM field components on the ID grid. 



to assign X y (i,t) and Y z (j, t) to the odd, respectively, even numbered lattice site, as shown in Fig. [I] for a grid of n 
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points. The equations ( |l5| ) and (|l6|) can now be combined into one equation of the form Eq. (JsJ) by introducing the 
n-dimensional vector 'Jj(i) with elements 

*(i,*) = |^ ) = ^?. ( y ) ' '° dd ■ (17) 
y ' \Y z {i,t) = ^E z (i,t), i even v ; 

The vector \&(i) describes both the magnetic and the electric field on the lattice points i = 1, . . . ,n and the ith element 
of ^f(t) is given by the inner product ^(i, t) = ef ■ *&(t), where e, denotes the ith unit vector in the n-dimensional 
vector space. Using this notation, it is easy to show that 

*(t) = E/(t)*(0) with U(t) = exp(tff), (18) 

where the matrix H is represented by two parts, 

H = Hi + H 2 , (19) 

with 

ri-2 

Hi = J2' h+U (e* ef+i - e i+1 ef) , (20) 

i=l 
n-2 

#2 = J! A+M+2 ( e 1+ l e f+2 - e l+ 2 e f+l) • ( 21 ) 

i=l 

Here, ftj = 1/ (5^/sijIj) and the prime indicates that the sum is over odd integers only. For n odd we have 

— *(M)=/? 2 ,i*(2,t) and — *(n,t) = -A l _i, n *(n-l,t), (22) 

such that the electric field vanishes at the boundaries (Y z (0,t) = Y z (n + l,t) = 0), as required by the boundary 
conditions Eqs. (||). 

The representation of H as the sum of Hi and i?2 divides the lattice into odd and even numbered cells. Most 
important, however, both Hi and Hi are skew-symmetric block-diagonal matrices, containing one lxl matrix and 
(n — l)/2 real 2x2 skew-symmetric matrices. Therefore, according to the general theory outlined in Sec. II, this 
decomposition of H is suitable to construct an orthogonal approximation 

Ui(t) = e THl e Tli2 (23) 

that is identical to the time-evolution operator U (r) up to first order in r. As the matrix exponential of a block- 
diagonal matrix is equal to the block-diagonal matrix of the matrix exponentials of the individual blocks, the numerical 
calculation of e rHl (or e rH2 ) reduces to the calculation of (n — l)/2 matrix exponentials of 2 x 2 matrices. The matrix 
exponential of a typical 2x2 matrix appearing in e rHl or e rH2 is simply given by 



exp 



a 



1 

-1 



\ / cos a sin a \ / *(i,t) 
*0'>*) J ~ I -sin a cos a J \V(j,t) 



(24) 



and represents the rotation of two elements of the vector *&(t) leaving all the other elements unchanged. This property 
of the time-evolution operator Eq. ( p3| ) provides the intrinsic possibility to parallelize the algorithms. Furthermore, 
it is even possible to alter the ordering of the products in the time-evolution operator U n (r) in order to construct an 
efficient implementation for a particular system. The plane rotations Eq. ( p4[ ) are performed by simply processing an 
arbitrarily ordered list S of pairs of EM field vector elements using 

Ui(t) = JJe 1 "*-^ * ^- e i e?) , (25) 
s 

instead of the odd-even decomposition Eq. ( p3| ) for which S = {(1, 2), (3,4), . . . , (n — 2, n — 1), (2, 3), (4, 5), . . . , (n — 
l,n)}. 

The implementation for ID can be readily extended to 2D and 3D systems, as has been illustrated in Ref. M. In 
2D, the TDME (|l]) separate again into two independent sets of equations and the discretization of continuum space 
is done by simply reusing the ID lattice introduced above. This is shown in Fig. || for the case of the 2D TM-modes. 
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FIG. 2: Positions of the three TM-mode EM field components on the 2D grid for n x — 9 and n y — 5. 




FIG. 3: Positions of the EM field components on the 3D Yee grid. 



The construction automatically takes care of the boundary conditions if n x and n y are odd and yields a real skew- 
symmetric matrix H. Correspondingly, in 3D the spatial coordinates are discretized by adopting the standard Yee 
grid jlfil , which also automatically satisfies the boundary conditions Eqs. (|§|). A unit cell of the Yee grid is shown in 
Fig. (f 

In general, the time step r and the distance S between next-nearest neighbor grid points are related due to the 
error that is introduced when the exact time-evolution operator U(t) is replaced by U„(t). We have ^ |[ 

\\U(r)-U n (r)\\ < j(d)(^^Y +1 . (26) 



Here, j(d) depends on the particular spatial discretization scheme used and a(n) represents the largest positive 
constant that appears as a prefactor in the exponential of the approximation C/„(r). We find a(2) — 1/2 from 
Eq. (10) and inspection of Eq. ( pd| ) yields a(4) = (1/2) (4a — 1) « 0.33. It follows that for a required spatial resolution, 



which determines the smallness of 5, the time step has to be chosen such that 



T < r* ee — - (27) 
a(n) 



in order to keep the error Eq. (26) small. As an example we consider a wave packet in a 2D cavity that is simulated 
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by a T4S2 algorithm. For numerical purposes we use dimensionless variables throughout this paper, where the unit 
of length is denoted by A and the vacuum light velocity c is taken as the unit of velocity, while the permittivity e and 
permeability \x are measured in units of their corresponding values in vacuum, respectively, Eq and hq. The cavity 
with e = l and /i = 1 is of size 19 x 15 and contains a dielectric medium with e = 2.25 and \x = 1 that has an inclined 
boundary. We plot in Fig. ^ the results of simulations in which the wave packet scatters on the dielectric medium. 
In the four pictures we show the EM energy density distributions that are obtained after simulation time t — 12.8 
for a fixed mesh size (5 = 0.1 and for four different time steps r. It follows from Eq. (p7[) that the upper limit for 



T=0.1 T=0.2 




FIG. 4: Energy density distributions at simulation time t = 12.8 for various time steps r obtained by the T4S2 algorithm for 
a fixed mesh size 8 — 0.1. The wave packet has been created for parameters (a x ,a y ) — (2, 1.73), (xo,yo) = (5,7.5), and k = 8 
(see for details Eq. (B4) in Sec. VB) and impinges on the dielectric structure from the left. 



the time step is given by t* = 0.3 in this case. For r = 0.4 the EM energy density distribution is, in fact, seen to 
change dramatically such that the results become meaningless. It should be noted that the limitation Eq. (E7|) on 
the time step is different from the Courant number which relates the time step t to the stability of finite-difference 
time-domain (FDTD) algorithms |l7J that are based on the Yee algorithm |u|. The algorithms presented in this 
paper are unconditionally stable by construction for any time step r and produce reasonable numerical results up to 
t = t* , a time step at which the Yee-based FDTD algorithms may have become unstable. 

We conclude this section by noting that our algorithms conserve the divergence of the EM fields only in ID systems 
but not in 2D and 3D systems. Although the initial state *&(t = 0) can be chosen such that the EM fields satify 
Eqs. (Q), the time- integration of the TDME by an algorithm based on the approximation U u (t) yields EM fields whose 
divergence quickly acquires a finite value and then remains constant in time. This is shown in Fig. |^ where we plot 
the computed norm of the magnetic field divergence in a 3D system as a function of time. The 3D system is an empty 
cavity (e = 1 and fj, = 1) of size 1.5 x 1.5 x 1.5 and we use the T2S2 algorithm. Though the divergence of the EM 
fields is not conserved in 2D and 3D systems, this error can be reduced by using smaller time steps or algorithms with 
higher-order time accuracy. This can be seen in Fig. |^, where we compare the algorithms T2S2 and T4S2 as a function 
of the time step r to show that the error in the EM field divergence vanishes for the TnS2 algorithm proportional to 
(rc/S) n . 



IV. VARIABLE GRID IMPLEMENTATION 



The poor man's implementation does not provide an optimal discretization scheme for physical systems of unregular 
geometrical shapes or with strongly varying permeability and/or permittivity. In a practical implementation of such 
systems the grid has to be variable with a small mesh size in one region of the system and a large mesh size in another 
region of the system. In this section we show how to implement a variable grid in such a way that the algorithms to 
solve the TDME remain unconditionally stable by construction. 

For the sake of simplicity we consider a ID system that is discretized using a variable grid as shown in Fig. In 
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FIG. 5: The norm of the divergence of the magnetic field in a 3D empty cavity (s = 1 and n — 1) of size 1.5 x 1.5 x 1.5 as a 
function of time t. The computation is performed with the T2S2 algorithm keeping the mesh size 5 = 0.1 fixed. 




0.1 



10 



5/cx 



100 



FIG. 6: The norm of the divergence of the magnetic field in a 3D empty cavity (s = 1 and n — 1) of size 1.5 x 1.5 x 1.5 as a 
function of 5/ct. The computation is performed with the algorithms T2S2 and T4S2. 



a straightforward implemention of the variable grid we would replace the constant next-nearest neighbor distance 5 
m Eqs. (|l|) and (H) of the poor man's implementation by the corresponding variable distance. It is convenient to 
write this substitution in the form 



1 



Si-ii - 8, 



i+l,i+2 



2A» i+i 



where <5jj is the distance between grid points i and j (see Fig. 0) and 

1 



(28) 



(29) 



is the averaged next-nearest neighbor distance. It can be easily checked that an implementation of the variable grid 
that relies on the replacement Eq. would destroy the skew-symmetry property of the corresponding matrix H 
(see Eq. (|T9|)) . This is unphysical: The original form of the Maxwell equations do have this property. However, a 
variable grid implementation that does preserve the underlying symmetry of Maxwell's equations can be constructed 
for a sufficiently smooth, variable grid. In this case, the second term in the brackets of Eq. (28) may be neglected and 
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FIG. 7: Positions of the two TM-mode EM field components on the ID variable grid. 



the replacement 



A 



= A, 



may yield a resonable approximation of Eqs. (15) and ( |16| ) for the variable grid implementation: 

d 1 



— X v (i,i) = — I -jf^M 







1 



Xy(i + 2,t) 



Xy{i,t) 



9t y/Si+l V Ai_|_l j i + 2 v //ii + 2 

The corresponding matrix H is seen to be skew-symmetric, 



T T 
e i+l e i _j_ e i+l e i+2 



-'i+2 c i+l 



(30) 

(31) 
(32) 

(33) 



and may again be separated into an odd and even part of which the exponents can be easily calculated following the 
same steps as given above in the poor man's implementation. 

It is obvious that this variable grid implementation can, in principle, be applied in any spatial dimension d. However, 
it is in general not possible to predict how to choose a grid that yields the best approximation to the true spectrum 
of eigenmodes of any non-trivial d-dimensional system. We therefore studied the criteria for the choice of suitable 
variable grids in particular systems numerically and present the results for a ID and a 2D system in the remainder of 
this section. 

The ID system under consideration consists of a cavity of length L = 10 with a constant permeability fi = 1 and 
a varying permittivity e. The permittivity deviates from its vacuum value (e = 1) due to the presence of a dielectric 
medium with e = 3 that is located in the middle of the cavity and extends over a length 2, as shown in Fig. ||. 
As a reference system we use a poor man's implementation with constant next-nearest neighbor distance 8 = 0.025 




FIG. 8: The ID cavity with the dielectric structure (solid line) and the two implemented variable grids: Aj^+i = {0.1 <-» 
0.05 <-► 0.025} (dashed line) and A M+i = {0.05 <-> 0.025} (dotted line). 



9 



and calculate the eigenmodes tu n of the corresponding matrix H. For two different variable grids we calculate the 
corresponding eigenmodes uj n and the deviation T(uj n ,uj n ) — 1 — oj n /u) n relative to the eigenmodes of the reference 
system. The two variable grids have in common that the dielectric medium and the transitions between s = 1 and 
e = 3 at both its sides is embedded in a grid of constant next-nearest neighbor distance which equals that of the 
reference system (A^j+i = S = 0.025). Furthermore, at the left end and at the right end of the cavity the next-nearest 
neighbor distance is constant over a length 2.5 and equals, respectively, A^j+i — 0.1 and A^j+i = 0.05 in the two 
variable grids. The transitions in the variable grids between regions of constant next-nearest neighbor distance involve 
abrupt steps between 



A M+1 = 0.1 <-» A M+1 = 0.05 «-» Aj.i+1 = 0.025, 
where we kept the intermediate distance Aj^+i = 0.05 over eight grid points, and between 

A M+1 = 0.05 <-> A M+1 = 0.025, 



(34) 



(35) 



respectively. 

In Fig. g we plot r(w n , Cj n ) for the first 50 eigenmodes of both variable grids. The relative deviation is seen to 
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FIG. 9: Relative deviation T(u> n ,£)n) for two variable grids. 



increase with the number of the frequency mode. As high mode numbers represent high frequencies this observation 
simply reflects the general fact that the accuracy of the eigenmodes depends on the smallness of the mesh size 
(numerical dispersion). Clearly, this also explains why the relative deviation T(u> n , uj n ) increases upto 2% for the 
variable grid with Aj.j+i = {0.1 <-> 0.05 «-* 0.025}, while for the variable grid with A^j+i = {0.05 <-* 0.025} this 
deviation remains well below 0.5%. For the first few frequency modes, however, we observe an increase in T(u n ,uj n ). 
This behavior can be related to the error that is introduced in the variable grid implementation by applying the 
approximation Eq. ([5(j|) instead of the exact replacement Eq. (|28|). To check this statement we plot in Fig. [l(] the 
deviation T(p, n ,w n ) for the first 50 eigenmodes of the two variable grids relative to the eigenmodes D, n that belong to 
the variable grids of the exact implementation Eq. (|2^) . We see that the increase of the relative deviation for the first 
few eigenmodes is, in fact, related to the error which is made by replacing the exact substitution Eq. ( pq) with the 
skew-symmetry conserving approximation Eq. (|30|). This approximation leads to oscillations of T(Q n ,u> n ) (and also 
r(w, w n )) that vanish with increasing frequency mode number. From extended numerical studies (results not shown) 
we find that these variations depend on several factors, such as the size in the difference between the largest and 
smallest distance Aj^+i of the variable grid implementation and on how abrupt A^j+i changes with i. In practice, 
it will be necessary to check the robustness of numerical results obtained by a variable grid implementation against 
small changes in its parameters. Although this may sound as a serious disadvantage, the next example of a 2D system 
shows that for realistic applications it may be by far more efficient to perform several simulation runs with a variable 
grid implementation than to use the poor man's implementation. 

The 2D system we consider is given by the L-shaped cavity depicted in Fig. [ll]. In order to satisfy the conditions 
Eq. (||) at the boundaries, the EM fields change very strongly close to the sharp edge of the cavity. Large spatial 
changes of the EM fields require a small mesh size. However, for the overwhelming part of the cavity a small mesh size 
would cause a waste of resources (computer memory and CPU time). Therefore, this system can be more efficiently 
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FIG. 10: Relative deviation r(f2»,cDn) for two variable grids. 




FIG. 11: The L-shaped 2D cavity with a variable grid (schematically). 



simulated by a variable grid implementation with an increasing number of grid points near the edge. This is done by 
a uniform increase of the number of grid points along both the x- and the y-direction as is schematically drawn in 
Fig. HU Furthermore, instead of using the odd-even decomposition of the time-evolution operator (corresponding to 
Eq. ( |23| ) for the ID system) on a square grid that would contain grid points outside the L-shaped cavity, we perform 
the plane rotations by processing a list S of pairs of the EM field vector elements at the grid points that actually 
belong to the L-shaped cavity (corresponding to Eq. ( pB| ) for the ID system). 

In Table | we present the results of a numerical simulation for the eight lowest TM eigenmodes in the cavity. We 
used the T2S2 algorithm imposing a poor man's implementation with S = 0.003125 and a variable grid implementation 



mode n T2S2 GdfidL 

constant variable constant variable 





grid uj n 


grid u n 


grid w„ 


grid w„ 


1 


2.9989 


2.9913 


2.9999 


2.9992 


2 


3.9807 


3.9500 


3.9740 


3.9720 


3 


4.9164 


4.8857 


4.9156 


4.9102 


4 


5.4150 


5.3843 


5.4077 


5.4004 


5 


5.5837 


5.5453 


5.5791 


5.5710 


6 


6.0592 


6.0209 


6.0580 


6.0494 


7 


6.7649 


6.7265 


6.7511 


6.7377 


8 


6.8876 


6.8492 


6.8797 


6.8674 



TABLE I: The eight lowest TM eigenmodes of the L-shaped cavity (see Fig. [ll| 
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constant grid S 


W3 


T (in %) 


0.1 


4.571 


7.5 


0.05 


4.740 


3.7 


0.025 


4.832 


1.7 


0.0125 


4.878 


0.78 


0.00625 


4.901 


0.31 


0.003125 


4.916 





variable grid A 


0.1 -» 0.05 


4.717 


4.2 


0.1 -» 0.025 


4.801 


2.4 


0.1 -> 0.0125 


4.840 


1.6 


0.1 -» 0.00625 


4.878 


0.78 


0.05 -» 0.003125 


4.886 


0.61 



TABLE II: Error in third-lowest eigenmode of the L-shaped cavity (see Fig. |ll| ). 



with a mesh size ranging from A = 0.05 to A = 0.003125. Very similar to the procedure described above for the ID 
system, the mesh size is decreased by a factor 0.5 and then kept constant for several grid points to smoothen this 
transition before the mesh size is decreased further. Our results are in good agreement with those obtained by the 
program package GdfidL |l8] for the same 2D system (see Table ||). In Table |j] we show the location of the arbitrarily 
chosen third-lowest eigenmode o>3 for several constant and variable grid implementations of the T2S2 algorithm. In 
all simulations we set S/ct — 10, where in the case of a variable grid 6 is replaced by the smallest mesh size. The 
relative error T of the frequency C03 is measured with respect to the frequency C03 = 4.916 of the system with constant 
mesh size 6 — 0.003125. The numerical results obtained within the variable grid implementation are in excellent 
agreement with the results of the poor man's implementation and the program package GdfidL. The T2S2 algorithm 
with the poor man's implementation and S = 0.003125 consumes about 150 times more CPU time and 10 times more 
computer memory than the T2S2 algorithm with variable grid implementation and A = {0.05 — * 0.003125}. Clearly, 
these numbers justify additional simulation runs that are required to check the robustness of numerical results against 
small changes in the parameters of a variable grid implementation. 



V. IMPROVED SPATIAL DISCRETIZATION IMPLEMENTATION 



Both conditional FDTD algorithms and the unconditionally stable algorithms TnSm suffer from numerical disper- 



sion due to the discretization of continuum space on a grid with a finite mesh size 17 . Methods to reduce numerical 



dispersion are taking a grid with a smaller mesh size or employing more accurate finite-difference approximations to 
the spatial derivatives. The former obviously can be also used in the poor man's implementation of unconditionally 
stable algorithms, however, for several reasons it may be more desirable to implement higher-order accurate approxi- 
mations of the spatial derivatives. For example, if one is interested in global features of the distribution of a system's 
eigenmodes, i.e. if we want to determine all eigenvalues, a higher-order accurate spatial derivative implementation 
would be strongly preferred. The computation of a system's eigenmode spectrum is performed by calculating the 
Fourier transform of the inner product F(t) = (0)\^f (t)) jj], [l9| [2(J. Using independent random numbers to initial- 
ize the elements of ^(0), the full eigenmode spectrum is obtained by averaging this Fourier transform. Taking just a 
smaller mesh size for the grid in the poor man's implementation does not only reduce the numerical dispersion but 
also gives rise to more eigenmodes with high frequencies. In order to obtain the eigenmode spetrum with the same 
spectral resolution, the sampling of F(t) would have to be done over smaller time intervals involving the computation 
of more data points. It is thus desirable to implement, instead, higher-order accurate approximations of the spatial 
derivatives that make a moderate use of computer resources in terms of CPU time and computer memory possible. 

The procedure for the construction of higher-order approximations to spatial derivatives is standard In the 

present case, we apply this procedure keeping in mind that Maxwell's equations (|^) are skew-symmetry and that 
the electric and magnetic field components are defined at particular grid points. The grid of a d-dimensional system 
with a constant mesh size of distance 8/2 between neighboring grid points is shown in Figs. ([l])-(|]). Without loss 
of generality we consider a ID system, where &(i,t) — ^(iS/2,t) is the ith component of the EM field vector and 
denotes an electric field compoenent for i even and a magnetic field component for i odd (see Sec. Ill for details). 
Applying the second-order accurate central-difference scheme the spatial derivative of the EM field component ^>(i, t) 
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is given by 



a_ n . t) _ « (t + i,o-«(i-i,o _ () + C(A (36) 



where ^/' 3 ^(i, t) = 9 3 ^(i, t)/dx 3 . Similarly, using the third- nearest neighbor EM field points at distance 35/2, we have 

a_ n . t) _ «( i + M)-« (i -^o _ 9^ (s)( . t) + C(J4) (3J) 

A fourth-order accurate approximation of the spatial derivative d^f(i,t)/dx is now constructed in terms of a linear 
combination of Eqs. ( |36| ) and (|37]) which is chosen such that the terms proportional to \&( 3 )(z,i) vanish. We obtain: 



8_ ni t) _ 9 + _ 1 ^(, + 3.0-^-M) ) +C(J4) (38) 

In practice, it is straightforward to implement the improved spatial discretization, since we can use the implementation 
of the central-difference scheme for the two terms separately and then combine the results according to Eq. (|3^). 
The corresponding matrix H of the ID system (see Eq. ([l9])) changes from tridiagonal to five-diagonal, but most 
importantly it preserves its property of being skew-symmetric. It should be noted, however, that the fourth-order 
accurate spatial derivative introduces errors at the boundaries since the calculation of d^^^) / dx for i — l,2,n — 1, 
and n refer, respectively, to grid points i = — 2, — 1, rt + 1, and n + 2 that lie outside the cavity and are implicitly 
assumed to be zero. 

It is obvious that the fourth-order accurate approximation of the spatial derivatives can be similarly applied in 
systems of any spatial dimension d. In the remainder of this section we study the numerical dispersion and the 
temporal and spatial accuracy of the algorithms for various ID and 2D systems. 

A. Numerical Dispersion 

We illustrate the difference in the numerical dispersion between the poor man's implementation and the improved 
spatial discretization implementation by a comparison of the eigenmode spectra of a ID empty cavity (e = 1 and 
(j, = 1) of length L. In ID, the continuum wave equation for the EM fields [Q, 

Id 2 d 2 1 

?5?-^J* M = ' (39) 

is solved by the ansatz ^(x,t) cx cos(ujt — kx + (f>) (with a phase (f> to distinguish electrical and magnetic field 
components) yielding the linear dispersion relation between frequency ui and wave number k: u = c\k\. Focusing on 
the effect of the spatial derivatives on the numerical dispersion, we assume perfect time integration of the algorithms 
and impose periodic boundary conditions on the EM field components: ^ p (i,t) oc cos(aj p t — k p S/2 + (f) with wave 
number k p = 2irp/L and —L/(25) < p < L/(2S). Applying the second-order accurate spatial derivative we obtain 

i) = 1 + 2, t) - 2* p ( i; t) + - 2, t)\ + 0{5 2 ) , (40) 
while for the fourth-order accurate spatial derivative we find 
d 2 , ,. s / 9 N " 



q^pW = \ T 5 ) [*„(i + 2,t)-2*p(M) + * p (i-2,t)] + 



2 



+ (^) [* p (t + 6,t)-2* p (* > t)+*p(i-6,t)] + 

+ 2, t) + *p(i - 2, t) - * p (i + 4, t) - * p (i - 4, *)] + 0(<5 4 ) . (41) 



245 / 
/ 9 



For m = 2 the analytical solution of the eigenmode spectrum for the mth-order accurate spatial derivative is given by 

to 2 = 2(^) 2 [l-cos(V)] , (42) 
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while for m — 4 we find 

2 3 

to/ = (|j ^ Q cos(Zfc p( 5) (43) 

1=0 

with coefficients C = 365/144, C x = -87/32, C 2 = 3/16, and C 3 = -1/288. We show in Fig. |l| that the dispersion 
relations, which we obtained numerically by the mth-order accurate spatial derivative implementation for a ID cavity 
of length L — 4, are in excellent agreement with the corresponding analytical solutions Eqs. (^) and (fl3|). It is 

35i ■ ■ ■ ■ ■ ■ 

CO 

30 



25 




FIG. 12: Numerical and analytical dispersion relations for the ID cavity of length L = 4 as obtained from calculations with 
mth-order accurate approximations of the spatial derivatives (m = 2,4). In both simulations we kept S = 0.1 and r = 0.01/c 
fixed. 

clearly visible that the dispersion relation computed by the poor man's implementation (T2S2 algorithm) suffers from 
numerical dispersion already at frequencies above u) — 10, whereas for a grid with the same mesh size the fourth-order 
accurate spatial derivative implementation (T2S4 algorithm) works well up to w = 15. 

B. Temporal and Spatial Accuracy 

To perform a systematic study of the accuracy of the algorithms as a function of the time step r and the mesh size 
S, we compute the difference between the normalized exact, and the approximate, \Jv m (i), EM field vector as 

obtained by the TnSm algorithm: 

A* n , m (t) = ||*(t) - *«, m (t)||. (44) 

We first consider the propagation of a Gaussian wave packet in a ID empty cavity (e = 1 and fi — 1) of length 
L = 30. At t = the Gaussian wave packet 

E z (x,t) = exp [-{x - x - ctf/a 2 ] (45) 

with standard deviation a = 2 is located at xq = 8. For t > the wave packet propagates with velocity c in the 
^-direction until it hits the right boundary of the cavity, becomes reflected, and propagates in the opposite direction. 
To derive an analytical expression of the exact EM field vector &(t), we expand E z (x, t) in the TM- modes 

E z [x,t) = —'^^a n s'm(nirx/L)sm(nir(xo + ct)/L), (46) 

n=l 

oc 

H y (x,t) = + a n cos(nnx/L) cos(nir(xo + ct)/L), (47) 

n=l 
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with coefficients 



a (rvK\- 



cxp ' -~ KT) 



(48) 



which ensure that the wave packet satisfies the boundary conditions Eq. (^) . Using Poisson's summation formula we 
find the following expressions for the EM field components 

E z {x,t) = ^2 [cxp(-(2nL + x + a;o-ci) 2 /CT 2 ) - cup (-(2nL + x - x + ct) 2 /a 2 )] , (49) 

71 — — OO 
OO 

H y (x,t) = ^2 [ ex P (~{2nL + x + x - ct) 2 /a 2 ) + exp (-(2nL + x - x + ct) 2 /a 2 )] , (50) 



from which the exact EM field vector is constructed according to Eq. ([Tt]) on the ID grid (see Fig. [I]). 

In Fig. [l3] we plot A\P„ )m (i) as a function of the simulation time t for fixed values of the mesh size 5 and the time 
step t using both the T4S2 and the T4S4 algorithm. We find that the error increases roughly proportional to the 




FIG. 13: The error A\E , n , m (£) as a function of the simulation time t for fixed values of the mesh size 8 — 0.1 and the time step 
r = 0.01/c. Results are shown for the T4S2 and T4S4 algorithm. 



simulation time: 

All/ n,m (t) = V2f n<m (r,6)t, (51) 

where we used the prefactor v2 to ensure that < f n ,m(T, 5) < 1. The linear dependence of A\I> Il!TO (t) on t is clearly 
visible only for the T4S2 algorithm but is also true for the T4S4 algorithm with a much smaller slope /4,4 (t, S) . Only 
at particular times t when the wave packet hits the boundaries of the cavity, the error A\&4 4 (t) is seen to increase 
nonlinearly in the time t and takes a value that is of the same order as A\&4 2 (£). This behavior, not described by 
Eq. (j5j), is present in fourth-order accurate spatial derivative implementations, in which the calculation of the EM 
field components close to system boundaries refer to several non-existing grid points. To study the error A 1 4'„. m (i) 
as a function of the time step r and the mesh size 5, we compute 

fn.m(T,5) = -i=i.A* n>m (t). (52) 

V2 at 

In Fig. |l| we plot /„, m (r, <5) as obtained for the ID cavity by the four algorithms T2S2, T4S2, T2S4, and T4S4 as a 
function of 6/ct for a fixed mesh size S. For each algorithm TnSm we find a linear decrease of log[/ njrn (r, S)] with 
increasing values log[(5/cr]. For the algorithms T4S2 and T4S4 we find that /4. m (r, 5) cx r 4 , while for the T2S2 and 
T2S4 algorithms f2,m( T ,3) °t t 2 - This numerical result is in agreement with the rigorous upper bound on the error 
of the EM field vector Eq. ([12]). For decreasing values of r, the error in the time integration becomes negligible small 
and f n .m(T,S) reaches minimum values which are indicated by the two lines 'Exact S2' for the algorithms TnS2 and 
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FIG. 14: f n ,m(T,S) as a function of 8/ct for the fixed mesh size 8 — 0.1. 
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FIG. 15: fn.mij, 8) as a function of 8 for the time step r = 0.15/c. 



'Exact S4' for the algorithms TnS4. In fact, these two lines represent the numerical results that are obtained for an 
exact time integration and mth-order accurate approximations to the spatial derivatives. 

Next, we study f n .m{T, S) as a function of the mesh size 5 for the time step r = 0.16 / c to ensure that the accuracy 
of the time integration remains constant. The numerical results are plotted in Fig. n5l We see that log[/„. m (r, 6)] 
decreases linearly with increasing log[l / 5] until it levels off. At this point, the total number of operations has become so 
large that it causes the numerical loss of accuracy. Outside this regime we find for the TnS4 algorithms f n ^(T, 5) cx 5 A 
and for the TnS2 algorithms f n .2(T,S) cx S 2 . In analogy to the upper bound Eq. (|l2|), the upper bound for the 
mth-order accurate approximation of the spatial derivatives is given by 



- *„, m (t)|| < Cn.mtS" 



(53) 



where C n ,m is a constant. 

We consider a second example to illustrate the numerical performance of the algorithms in 2D systems. For the 
initial wave packet in the 2D cavity we make the ansatz 



E z (x,y,t) = sin(k(x - x - ct))exp[-((x - x a - ct)/a x ) 10 - ((y ~ yo)/<J v ) 2 



(54) 



At t = the wave packet is centered at {xq, yo) and moves at t > with velocity c in the x-direction. The energy of the 
wave packet is fixed by the wave number k in the oscillating factor and its envelope is Gaussian along the y-direction 
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and has sharp edges along the x-axis (due to the exponent 10). 
contains two objects with dielectric constants e = 5 and [1 = 1. 
(&x,& y ) = (1.66,1.29), (xo,yo) — (3.5,5.5), and k = 5. In Fig. 
T2S2 and T2S4 algorithms with different mesh sizes relative to a reference EM field vector *f?(t) that is obtained from 
the T2S2 algorithm at mesh size S = 0.025. In all simulations we kept r = 0.15 /c fixed to compare measurements of 
constant accuracy in the time integration. In Fig. |l^ we show (a) the energy distribution of the initial wave packet 



The 2D cavity of size 12 x 10 with e = 1 and \i = 1 
The parameters of the propagating wave packet are 
16 we show the results for the error Eq. (H) of the 




e 



FIG. 16: The error AtD„, m (r,t) for various mesh sizes and algorithms in 2D. (a) Initial energy density distribution \P(r,t) 2 . (b) 
Reference energy density distribution ^2,2 (r, t) 2 at t = 6 using the T2S2 algorithm with 6 = 0.025. (c) The error Aa)2,2(r,t) 
on the energy density distribution at t = 6 using the T2S2 algorithm with 8 = 0.1. The relative deviation is 26%. (d) The 
error Aw2,2(r, t) on the energy density distribution at t — 6 using the T2S2 algorithm with S — 0.05. The relative deviation is 
5.9%. (e) The error Aw2,4(r, t) on the energy density distribution at t = 6 using the T2S4 algorithm with 8 = 0.1. The relative 
deviation is 5.1%. 



(t = 0) and (b) the reference energy density distribution after simulation time t = 6 using the T2S2 algorithm. In 
(c)-(e), the normalized spatial distribution of the error in the energy density distribution, 

A«w(r,i) = |*(r,t) 2 -*„, m (r,t) 2 |, (55) 

is shown for, respectively, the algorithm T2S2 with 5 = 0.1, the algorithm T2S2 with 5 = 0.05, and the algorithm 
T2S4 with 6 = 0.1. We find that the improved spatial discretization implementation T2S4 with 6 = 0.1 performs as 
well as a poor man's implementation T2S2 with half the mesh size. The main advantage of using the T2S4 algorithm 
is that it used only 20% of the computer memory and 10% of the CPU time with respect to the T2S2 algorithm. 



VI. CONCLUSIONS 



We have demonstrated that the previously introduced family of unconditionally stable algorithms to solve the time- 
dependent Maxwell equations can be implemented with a grid of variable mesh size and with a fourth-order accurate 
approximation to the spatial derivatives. The performance of the algorithms has been shown to increase significantly as 
compared to the previously applied poor man's implementation while at the same time their property of unconditional 
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stability by construction is preserved. Performing numerical simulations on various physical model systems, we 
found that a variable grid implementation can save orders of magnitude in computer memory and CPU time for a 
physical system of unregular geometrical shape or with strongly varying permeability and/or permittivity. Similar 
enhancements have been obtained for the fourth-order accurate spatial derivative implementation which does not only 
reduce the numerical dispersion but also improves the temporal and spatial accuracy of the algorithms significantly. 
Clearly, in close analogy to the implementation of the fourth-order approximation of the spatial derivatives, the 
algorithms may be improved by constructing higher-order approximations. In general, we conclude that the family of 
unconditionally stable algorithms does not only preserve the fundamental symmetries of the time-dependent Maxwell 
equations but is also characterized by a high degree of flexibility that allows to construct implementations that are 
required in different kinds of specific applications. 
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